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■ We analytically calculate the equilibrium sequence of the corotating binary 



stars of incompressible fluid in the first post-Newtonian(PN) approximation 
of general relativity. By calculating the total energy and total angular mo- 
mentum of the system as a function of the orbital separation, we investigate 



^ ■ the innermost stable circular orbit for corotating binary(we call it ISCCO). 



It is found that by the first PN effect, the orbital separation of the binary 

at the ISCCO becomes small with increase of the compactness of each star, 

and as a result, the orbital angular velocity at the ISCCO increases. These 

behaviors agree with previous numerical works. 
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I. INTRODUCTION 



The laser interferometers such as LIGO §, VIRGO 0, GEO600 § and TAMA300 § 
are currently being constructed and the detection of gravitational waves is expected in this 
decade. One of the most important astrophysical sources of gravitational waves for these 
detectors is the coalescing binary neutron stars(BNS's) because the gravitational waves 
emitted in the inspiral phase (the so-called last three minutes) |j] have frequencies in the 
sensitive region of these detectors, i.e., from 10Hz to 1000Hz. We will be able to know each 
mass, spin, and so on of the BNS's if we could obtain an accurate theoretical template for 
data analysis ||. Hence, much theoretical works have been done to complete it f| ||. 

When the orbital separation of the BNS's becomes a few times of the neutron star(NS) 
radius as a result of the radiation reaction of gravitational wave emission, the hydrodynam- 
ical effect becomes important. In such a phase, the wave form of gravitational waves is 
expected to be sensitive to the NS structure, especially the relation between the radius and 
mass of the NS. Therefore, if gravitational waves from such a phase are detected, we may 
constrain the equation of state(EOS) of the NS JK| [|Il|] |12|]. In particular, the important 
quantity is the location of the innermost stable circular orbit (ISCO), which will have an 
information on the EOS of the NS. 

There have been many analyses for the ISCO using some approximations of general rel- 
ativity [|13|] [|14| |L5| fl6|], but in order to determine the precise location of the ISCO, we 
need fully general relativistic(GR) numerical simulation, which is a very difficult method. 



Recently, Wilson and his collaborators |L7j have performed numerical simulation for obtain- 
ing the equilibrium sequence of the BNS and its ISCO solving semi-relativistic equations. 
The results they have obtained are very interesting, but the accuracy of their results is still 
in question because they do not show any calibration of their numerical code using test 
problems. Furthermore, they do not seem to perform any important analysis to their nu- 
merical results. When we carry out a large numerical simulation, it is required to perform a 
detailed analysis after the computing in order to explain the numerical results. The analysis 
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is desired to be done by comparing the analytical or semi-analytical estimates. If the nu- 
merical results qualitatively agree with such an analytical calculation, we can firmly believe 
the numerical results and also can understand details. Hence, with a large scale simulation, 
it is favorable to prepare some analytical models in order to understand the essence included 
in the numerical results. 

For that purpose, we here analytically solve the Darwin problem in the first PN ap- 
proximation. The Darwin problem is concerned with the equilibrium and the stability of a 
homogeneous fluid star rotating around another one taking into account the mutual tidal 
interactions fl8 |. We assume that each NS in the binary system is composed of the in- 



compressible and homogeneous fluid. The reward is that all the calculations can be done 
analytically. This means that we can obtain a strict solution of a BNS including the GR 
effect without any large supercomputing. Our results will be very helpful for understanding 
how the finite-size effects as well as the GR one of the NS's influence to the location and 
the orbital angular velocity at the ISCO, and for the analysis of large scale simulation. 

This paper is organized as follows. In section II, we show the basic equations to solve 
the PN Darwin problem. By using the first tensor virial(TV) equations, we derive the 
angular velocity of corotating binary systems in section III. In section IV, we show the 
equations of the total energy and the total angular momentum for corotating binary systems. 
In section V, we calculate the equilibrium sequences and determine the location of the 
energy and angular momentum minimums(we call it the innermost stable corotating circular 
orbit (ISCCO) to distinguish it from the innermost stable circular orbit (ISCO)0). Section VI 
is devoted to summary. 

Throughout this paper, we use the unit of G = 1, and c denotes the light velocity. Latin 
indices i, j, k, ■ ■ ■ take 1 to 3, and Sij denotes the Kronecker's delta. We use Iy and 4-^ as 



^This point is the secular instability limit, and not the dynamical instability limit(i.e., ISCO). 
The ISCO will locate inside the ISCCO E3|. 
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the quadrupole moment and its trace free part, 



1 3 

6 k=l 



of each star of binary. 



II. FORMULATION 

Non-axisymmetric equilibrium configurations of uniformly rotating incompressible fluid 
in the first PN approximation are obtained by solving the integrated form of the Euler 
equation and the Poisson equations for gravitational potentials consistently. The integrated 



form of the Euler equation was derived by Chandrasekhar and it can be written as [20 



21 



P 1 (P\ 2 



p 2c 2 \ p 

C 2 I 2 ' C 2 V ' rvj I " ' 4c 2' 



U ~ + i + \ ( 2 ^ U ~ X n + K) ) ^ + TT^ + const., (2.1) 



where 



w2= + (2 ' 2) 

and U, X , Xq, and (3 V are the gravitational potentials. In this paper, we consider the 
equilibrium sequences of BNS's of equal masses(M 1 = M 2 = M) whose coordinate separation 
is Rf^ We assume that the center of mass of a star (star 1) locates at the origin and the 
other one (star 2) locates at (xi, x 2 , £3) = (— R, 0, 0)(see fig. 1). Due to the symmetry, we 
only pay attention to the equilibrium configuration of star 1 in the following. 

The main purpose of this paper is to calculate the PN correction of the angular velocity, 
the energy and the angular momentum for the Darwin problem. In the Newtonian order, 
the angular velocity(f2N) becomes ^3[(see below for derivation) 



=The coordinate condition in this paper is the standard PN one |22] 
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Thus, in the PN approximation, we can expect that the following types of quantities will be 
the main terms in the PN order: 

MM MM Mai M , Ma^ M , . 

~ — x -, ~ — x — -, ~ — x -, and ~ — x — -, 2.4 

R 3 a c 2 ' R 3 Re 2 ' R 5 a c 2 ' R 5 Re 2 ' K 1 

where ao is a typical radius of the star, and we use the relation 4-n ~ Ma\. We will derive 
four types of terms shown above and the correction of the energy and the angular momentum 
by them below. 

Since we consider the incompressible fluid, the gravitational potentials inside each star 
are expressed as the polynomial form of the coordinates Xj. For the purpose of obtaining 
the PN corrections shown above, we need to take into account the coefficients of the terms 
such as x^x^x™ 3 in X , Q 2 Xq, ^l 2 (3 ipi and so on, where < m 1 ,m 2 ,m 3 < 5 and < 
mi + ni2 + ffl3(= Trit) < 5, up to 0(R~ 5 ) for the case m t is odd and up to 0(R~ 3 ) for 
the case m t is even. In the following subsections, we solve the Poisson equations for the 
gravitational potentials to derive such terms. 



A. Newtonian Quantities 

Each gravitational potential is composed of two parts; one is the contribution from star 
1 and the other is from star 2. In the following, we denote the former part as 1_>1 , and the 
latter one as 2- * 1 , where denotes one of the potentials. We also define (p 1 ^ 2 and 2 ^ 2 as 
the contribution from star 1 to 2 and star 2 itself, respectively. 



Following previous authors |L3| ||23|| , the configuration of each star of binary in the 
Newtonian order is assumed to be an ellipsoidal figure of its axial length a%, a 2 and 03. In 
this case, the solution of the Poisson equation for the Newtonian potential 

AU = -47rp, (2.5) 

is written as U — U 1 ^ 1 + U 2 ^ 1 , where 



and 



M j xi 2x\ ~ x l ~ x i 
~R 1 ~~R + 2R 2 + 



(2.6) 



8xf + 3z§ + 3x| - 24x 2 (:r 2 , + + 6x|x^ ' 



8P 4 



3^11 / 3xi 12xf — 5x. 
+ 1 2B Y [ + 2~R2 



2 _ r 2 X 

2 dX 3 



+ 



2i? 5 



(^22^2 + ^33^3) 



(2.7) 



(2.8) 



Aij... are index symbols introduced by Chandrasekhar [18], and A = J2iAiaf is calculated 
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A = aia 2 a 3 / 
Jo 



du 



yj \a\ + u)(al + u)(a§ + u) 



a? A), 



(2.9) 
(2.10) 



where 02 = 02/^1 and 0:3 = a^/a\. Note that U 2 ^ 2 and [7 1- * 2 are obtained by changing x\ 
into —{x\ + i?) in t/ 1 ^ 1 and U 2 ^ 1 , respectively. 

In the Newtonian order, the pressure is written as 



(2.11) 



Pn is calculated from the scalar virial relation as 



Po 



npAo - (a? + 



M 
2~R3 



o 2 2 2 

zdj — a 2 — a 3 



+ 0(R- 5 ), 



(2.12) 



and «2 and «3 are determined from 



pa 2 

pal 
fo_ 
pal 



n 2 N m 

^ 2 2P 3 



M 
2~R3' 



(2.13) 
(2.14) 
(2.15) 
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B. X, 



o 



As in the case of U, the PN potentials are divided into two parts as Xq = X^ 1 + X^ 1 , 
and we consider them separately. 

• Contribution from star 1: 



Xq is derived from the Poisson equation fl20 



AX*-* 1 = 4vrp 



2vrp (a - Mx] 



3P n 



and the solution becomes 



2U 



2->l 



(2.16) 



Xq 1 ^ 1 = -aoU 1 ^ 1 + a 1 D 1 + 5>,A 



-^(2D n - D 22 - D 33 ) - ^(-2D U1 + 3,0x22 + 3 Ass 



M 



R? 



i? 4 



(2-17) 



where 



3P n 2M 34- 



a = 2npA H h n 

p it 

2M 9Pn 
pa/ 



A, Ai, and A« are the solutions of equations 



ii 



R 3 



(2.18) 
(2.19) 
(2.20) 



AA = 


-47rpxi, 


(2.21) 


AA< = 


-47rpx-, 


(2.22) 


AA« = 


— 47rpxix^ , 


(2.23) 



and the solutions at star 1 are |T3 



[2.24) 



Dm = up 



+ ~A ( B i - 2 E B ^ 2 l + E E Bum^xA 
4 \ j i m J 

i \ ^111 -E^ni^i) x ? 

+ -a{ (B u - 2 £ B m x 2 + £ £ B nim x 2 x 2 m ) 
Dm = irp a\a* ^A U i - E ^uu^ij x x x 2 

+ \ a i a2 i [ B li ~ 2 E ^Mi^ + E E B lilmXfx 2 m \ X 1 



where B^... are index symbols defined by Chandrasekhar [T8 



Contribution from star 2: 



The equation for Xq is 

AX 2 ^ = Anp ^p (A - £ ^y?) + ^ (l ~ E f ) + a*/ 1 " 
where ?/i = — (xi + -R), y<i = x 2 and ?/ 3 = x 3 . The solution is written as 



2^1 



M 



where Dfj?} are calculated from the same equations as the case of D^... [|T^]. 

AD 2 - 1 = -ATxpy iyj 



The solutions of Df* are 



Di 



2->l 



2-^1 



where 



hi ( 1 2xi | Gxf 3^2 3x| ^ o^r- 3 ^ 



R 2 \ R 

h2 fX 2 _ 3X1X2 

B?\~R~ R 2 



2R 2 

+ 0(R~ 3 ] 



-2x? + 3xi(x2 + x 



R 



R 



2R 2 



2R 3 



+ 0(R- 



— hill — o E 



s 



C. X n 



Contribution from star 1: 
The equation for X^ 1 is 



inp \x\ + x\ + Rxi + 



R 2 



(2.35) 



Then the solution is 



X, 



l-t-i 



R2 Trt-*l 



-2\D 11 + D 22 + RD l + — U 



Contribution from star 2: 

The equation Xq^ 1 may be written as 



lirp (yj + 2/2 - Ryi + ^- 



Then, using P 2 p.\ the solution is easily derived as 



X, 



2->l 



-2 P 



ii 



2 



-^22 


- PP 2 ~ 


1 + ?^ 
4 


2/n 


2/22 




P 2 Xl 


P 


(-1) 



(2.36) 



(2.37) 



(2.38) 
(2.39) 



D. p„ 



The definition of /3„ is E 



7 / P 
= -- (xiPi + x 2 P 2 + -Pi 



xi + - J P 2 ,2 + a%Pi,i - \ xi + -J x 2 (P 1>2 + Pj 
where Pi and P 2 satisfy, 

APx = -47Tp + D , 

-Anpx 2 . 



AP 9 



(2.40) 



(2.41) 
(2.42) 
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Pi is also written as P.- 1_>1 + P; 2_>1 , where 



and 



Pl- l = D l + ^U^\ (2.43) 
P^ 1 = D 2 , (2.44) 



2 

P 2 2 ^ = D\^\ (2.46) 



E. The Collection 



Substituting expressions for the gravitational potentials derived in subsections A— D into 
Eq. (p.l|) , the integrated form of the Euler equation is written as 

— = u + su + -w 2 n 2 

P 2 



1 



2 2 
T 



l>m 



l>m 



+const. 



(2.47) 



where <5P is the PN correction of U which we will mention in the next section. In the 
following, we do not need 70, 7«, and 7^ (see below), but need 0q, 0i and which are 



0o 



Pi 



Mirp ( 6P 1L4 3 2 , 2 , \ 9M 2 
R 2 { 5vrp 2 10 2 1 2 / 4P 3 



+ 



P 4 



1 9 / 3Pn\ 

: 7rp^n(-aiAi - 2a\A 2 + A ) + - J2[ 2 *P A i + — *uu 

1 1 \ P a i ) 

3 9 / 3P 

+ -7rpMa 2 1 (-2a 2 1 B 11 + a 2 P 12 + a 2 P 13 ) - -n P P u I 2A + — | 



+— (118/n - 93/ 22 - 59/33 , 
8P 5 

Mvrp /3 2 . 2 . 1 , 



+ 



P 4 



M 



(2.48) 
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02 



03 



011 
022 
033 
012 
013 
023 



Mnp 

1 

+ 1J 



--irpMal(— 4a^m + 6a 2 £> m — 3a\Bu 2 — "&a\B 

1 3 
-a\A l2 - a 2 2 A l2 + 3a|A 2 2 + A 1 - - 



M 2 



Mvrp 
1 



M?rp 
Mvrp 
Mvrp 
Mvrp 
Mvrp 
Mttp 
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-7rp^ii(-a 2 Ai2 - 2a 2 ,/4i 2 + 6a 2 2 A 22 - 3A 2 + 2A{) 
3 

--irpMa{{-2a{B n2 - 2a 2 A 122 + a\B l22 + a^ 12:> , 

3 1 N 

-a\A l3 + a 2 2 A 23 - -A 3 



+ 



29M 2 
8R 5 ' 



(-a 2 /4 13 + a 2 /4 23 - -A 3 ) 

^7rp^ n (a 2 /4 13 + 2a 2 /4 23 - A3) - + \ \27rpA + 

3 1 39M 2 

--TrpMa 2 (-2a 2 B 113 - 2a^ 133 + a\B l23 + a 2 5 133 ) + 



M 



3 

--7rpMa 2 (-2a 2 5 113 - 2a^ 133 + a 2 2 B 123 + a 2 3 B l33 ) 
2a\Ami + -a 2 (— 2a 2 £? mi + a2-B 1112 + a 3 -B m3 



4 

3 

1222 + ^i( — 2a 1 -Bn 22 + a 2 B 1222 + a 3 B 



' 3 D\ 223l 



—3a 2 a 2 A 

—3alalA 1333 + -a 2 (— 2a 2 £? 1133 + a2-B 1233 + a 2 £?i 333 ) 



3 

2a\A lll2 - 3a\a\A n22 + -a 2 (-2a 2 5 1 



1-Blll2 + a 2 -Bll22 + a 3 -Bll23) 
3 

2a?A m3 - 3a 2 a 3 A n33 + -a 2 (-2a 2 J B m3 + 03.81123 + a 2 3 B U33 ) 

3 

-3a 2 (a2v4i 223 + a 3 Ai 233 ) + -a 2 (-2a 2 Si 123 + alB 1223 + a 2 3 B 



'1233) 



(2.49) 



(2.50) 



(2.51) 
(2.52) 
(2.53) 
(2.54) 
(2.55) 
(2.56) 
(2.57) 



III. THE POST-NEWTONIAN ANGULAR VELOCITY 



In this section, the orbital angular velocity in the PN order is derived using the first TV 
equation. The first TV relation is derived from 



d x 



dP 

dxi 



0. 



(3.1) 



Substituting Eq. fl2.47p to Eq. (|3.1|), we have Eq.(|2.3|) for the Newtonian order. In the PN 
order, the explicit form becomes 
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= ^M5n 2 + J d 3 xp^n 2 N + M(3 + + (3 2 I 22 + (3 3 I 33 + 5/3 n / ull 

r dU 

+3/3i 2 Al22 + 3/313/H33 + /922/2222 + ^33^3333 + /?23^2233 + ^ j d?Xp—, (3.2) 

where £1 is the x\ component of the Lagrangian displacement and £f2 2 denotes the PN 
correction of the orbital angular velocity. 

A. The Post-Newtonian Deformation 

The density profile of an incompressible, homogeneous sphere in the PN order is the 
same that that in the Newtonian order. However, the density profile of an ellipsoid in the 
PN order is different from that in the Newtonian order. The Chandrasekhar's method to 



calculate the PN effects on the uniformly rotating isolated bodies is as follows [jlJJ. First, 
he constructs the ellipsoidal figures in the Newtonian order. Then, he regards the PN effect 
as a small perturbation to the Newtonian configuration, and calculates the deformation 
from the Newtonian ellipsoid using the Lagrangian displacement. Solving the equations for 
the Lagrangian displacement and calculating the correction to U by the deformation, he has 
obtained the equilibrium configuration in the PN order. In this paper, we follow his method. 

Although the deformation from the ellipsoidal figure will really occur by the PN effect, 
we will show in this subsection that it is a small effect compared with leading terms shown 
in Eq.(2.4). Hence, for readers who are not interested in such a verification, we recommend 
to skip to the next subsection. 

In order to express Eq. (|2.47| ) using only the coordinates and the index symbols, we have 
to calculate SU which is defined as 

5C/ = - V A / d 3 x' T -^—. (3.3) 

Here, & denotes the Lagrangian displacement of the fluid element induced by the PN cor- 
rection, and is used to guarantee that the pressure of the deformed Darwin ellipsoid at the 
boundary becomes zero. This means that it is determined from 70, ji, 7^, fto, Pi and /% in 
Eq. fl2.47D in general. 
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In choosing £ i; we require it to be divergent free because 5p = — J2i °i(pii) ■ To obtain 
5Q 2 up to 0(R~ 3 ) x Q^, we only need the following^ 

d 0) =^,0,0V (3.4) 



2 



t(21) 



= Qx^O, -arjars^ , (3.5) 
• 22) = (0,xix 2 , -xix 3 ). (3.6) 



Thus, we set 



e t = s d 0) + s 21 d 21) + «f 2) , (3.7) 

where So, £21 and S22 are constants. In this case, we can get SU of the contribution from 
star 1 as 

d£7 = vrpAiSbzi - -— — S 2 i - -5 5 — ^22, (3.8) 

\2 ox 1 0x3 j \ 0x2 ox 3 J 

where 

D u = npa\a\ ^A u -J^Amx^j x x x { (i ^ 1). (3.9) 
The contribution from star 2 is defined as 

^^-T.i-Ux'^-. (3.10) 



dx{ 

Here, the integral is performed in star 2, and £| is defined for star 2. Note that is the 
same as |^|, but £• has opposite sign to & if components of & are the even function of X{. 
Using this definition, we have 



**In ref. [19| where Chandrasekhar calculated the PN configuration of the Jacobi ellipsoid, he 
introduced other type of displacements, such as £, k (k = 1 — 5). Coefficients of these displacements 
are, however, of 0(i?~ 3 ) in this paper because the ellipsoidal displacements of a star are generated 
by the small spin of each star and by the tidal force of the companion star both of which deform 
the stars only by 0(R~ 3 ). Then, the displacements will contribute to the orbital angular velocity 
from a higher order of R^ 1 as x 0(i?~ 5 ). This is the reason why we do not need 70, 7« and 7«. 
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1 BU 2 ^ 1 fldD 2 ^ 1 8D 2 ^ 1 \ fdD 2 ^ 1 BD 2 ^ 1 \ 

5U 2 ^ = + l 1 -^- - ) S 21 + l^- - ^2- ) S 22 , (3.11) 

2 axi \2 axi ax 3 J \ ox 2 ox 3 J 

where 



D '< =-^U--^J (312) 

In order to get the orbital angular velocity, we also have to calculate the second and the 
last terms of the right-hand side of Eq. ( |3.2|) which are occurred by the displacement of the 
fluid element by the PN correction. The second term is 

y^(s + S 21 |V (3.13) 

and the last one can be evaluated as follows: First, the contribution from star 1 is zero 
because 

OU 1 1 f 3 / \ f 73 ///\ %i — x[ 



r/|3 



5 / d 3 xp— = -5 i d 3 xp(-x) [ cPx'pOx.' . 

J OX\ J J |x — x' 

= - / d 3 xp(x)5 f rfVp(x') ^ ~ X '\ + f d 3 x'p(x')5 f d 3 xp(x) f 1 - X \ 
J J |x — x | d J J |x — x/| d 

= 0. (3.14) 
The contribution from star 2 becomes 

The first term on the right hand side of Eq. fl3.15|) denotes the force which the displaced 
element of star 1 receives from the non-displaced potential of star 2. On the other hand, 
the second term denotes the force which the non-displaced element of star 1 receives from 
the displaced potential of star 2. For the Lagrangian displacement given in Eqs.(|3.4|) — (|37^), 
the first and the second terms of Eq. ( |3.15| ) are equal, and using Eq. ( |3.11| ) , we can calculate 
Eq. flgnSD as 



&xp% = 2 / d 3 xp^- , (3.16) 

oxi J ax\ 

2M ( 18^ii \ , _ 2M / , 9^ mi , 9^ u /n , 12J 1133 
— bo-=-r- M H — — + in 



R 3 V R 2 J R 3 V R 2 MR 2 R? 

2AM 

+'5*22 n c (-^1133 — -^1122)- (3-17) 

R* 
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At the deformed surface of the Darwin ellipsoid, the pressure must vanish. This condition 
becomes 

(-) =-2*Y,^r + tU + W*<* 
\pJs P i a i 2 

+70 + lltf + T/m^m + X 1 A) + Y Pl X f + H Plmtf x li + COnst. 
I l>m \ I l>m J 

= 0, (3.18) 

where we use the equation of the boundary surface at the deformed Darwin ellipsoid 

S = E4-l-2£^ = 0- (3.19) 
l i i i 

Eq.( |3.18|) must be satisfied at the original surface of the Darwin ellipsoid, ^xf/a 2 = 1- 



ower 

122; 



Coefficients, So, S21 and S 22 , are determined from the fact that the terms of the odd p 
of xi must vanish. The coefficients of x±, xf, X\x\ and X\x\ in Eq.( |3.18| ) (Qi, Qui, Q 
and Q133) are as follows: 

Qi = (~— 2 + *pAi + 5s )so + + ft 

\ paf 2 J 2 

+S 21 Ttpa\ (-alAu + ~B n + a 2 3 A 13 + 0(iT 3 )) - S 22 ixpa\ (a\A l2 - a 2 3 A 13 ) , (3.20) 

Qui =pi + S 21 irpal ( ^-7 + 2a\Am - a 2 3 A 113 - ^B ni ] 

V V«i 2 / 

+S 22 7ipal (a 2 2 A ll2 - agAns) , (3.21) 
Q122 = fa + S 21 ixpa\ (a\A 112 - ^ B U2 ~ a\A X2 z 



+S 22 npaj ( 2 f° + 3ajA 122 - ajA 123 ) , (3.22) 

Q133 = P3 + S 21 irpal ( 2P ° 2 2 + a\A ll3 - \b 113 - 2>a\A 133 ] 

\irp 2 aiai 2 ) 

+S 22 npa 2 1 ( 2 f + a\Ai 23 - 3a 2 3 A X33 J . (3.23) 
\irp 2 aiaj ' J 

The conditions to determine S , S 2 \ and S 22 are 

Ql + <5l33«3 = 0, (3.24) 

Qlllfll - (513303 = 0, (3.25) 

Ql22fl2 - Ql33«3 = 0. (3.26) 
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It is found that the coefficient of So in Q\ vanishes when we substitute Eg. (|3~2|) to Eq.( |3.24|) 
and remove 5fl 2 . Thus, Sq is indeterminant, and the first equation becomes trivial one. 
This is also verified by calculating 3Qma 2 + Qi 22 a 2 + Qi33 a h which coincides with —5Q±. 
This fact seems to be very natural because Sq is a simple transformation of the center of 
mass. As a result, S 2 i and S 22 are determined by the second and third equations. From the 
combination 2 x ( [3.25D — (|3.26|) , i.e., 2Qma? — Qi22«2 — Qi33 a 3 — 0, we obtain 

= 2fi x a\ - f3 2 a 2 2 - f3 3 a\ 



+S 2 iirpal 



4P f 



o 



., „ + 4a 1 yini - a x a 2 A 1V2 - 3a 1 a 3 A 1 i 3 + a 2 a 3 A 123 

+303^33 - a\B ni + -{a 2 2 B 112 + a\B n3 ) 
+S 22 npa 2 1 \2a\{a 2 2 A 112 - a 2 3 A 113 ) - 3ajA 122 + 3ajA 133 

where we omit higher order terms of R~ l . 2f3ia\ — j3 2 a 2 — (3 3 a 3 is calculated as 



Mnp 



R 2 

- o(i?- 4 ), 



5 15 
2a\A x - -a\A 2 + -a\A 3 + ~a\(a\A X2 - a 2 3 A 13 ) 



0{R 



-4\ 



(3.27) 



(3.28) 



where we use the relations 



a 2 2 A 2 = a\Ai + 0(R~ 
a jA 3 = a\A x + 0(R~ 
a\a\A X2 = a\A 3 + 0(R~ 
a\a 2 3 A 13 = a\A 2 + 0(R~ 



(3.29) 
(3.30) 
(3.31) 
(3.32) 



The coefficient of S22 is also of 0(R 3 ) because a 2 A XX2 — a 2 A XX3 = 0(R 3 ) and a 2 A X22 — 
a ±A l33 = 0(R- 3 ), but that of S 2 i is 0(1). 



On the other hand, from Eq. fl3.26|) , we have 

= f3 2 a\ - f3 3 a\ 
+S 21 irpa 2 1 



' a\{a\A XX2 - a 2 3 A X i 3 ) - -ji.a 2 2 Bu 2 - 03.6113) - a 2 2 a 2 3 Ai 23 + 303^133 



+S 22 irpa\ 



irp 2 a\ 
4P 



, + 3a 2 Ai22 - 2a 2 a 3 Ai2 3 + 3a 3 Ai 33 
7rp z a| 



(3.33) 
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The coefficients of S21 and S22 are of 0(1) in this equation, and we can also show f3\a\ — 
f3 2 a\ = 0(R~ 4 ). Thus, we can conclude that both S21 and S22 are of 0(R~ 4 ). Accordingly, 
we do not have to consider these Lagrangian displacements in the following. 

B. The Post-Newtonian Orbital Angular Velocity 

Since S is arbitrary, we set S = 0. Hence, Eq. (|3.2|) is calculated as 
R 



-M5Q Z = M(3 + 3/hhi + P2I22 + Pshs + + 3/? 12 / 1122 

+ 3/3i3/n33 + /?22-^2222 + /?33^3333 + $23-^2233 

M 2 A 9M 3 84 . M 3 . 9 , , 9x 

27tpA + - —jTrp^nA) + 7777^( 28a i ~ 14a 2 ~ 9a 2 ,), 



(3.34) 
(3.35) 



i? 2 r u 4i? 3 5i? 4 r 11 u 10i? 5 

where we make use of the relations A = J2i ^i a h Eqs. (|2.13|) — ( p.!5| ), equations for Iujj^\ and 
reduction formulas for A^... and B^... shown in the section 21 of Chandrasekhar's textbook 
I8| . Then, Q 2 is written as 



n 2 



2M 



R 3 
18^ii 



If , 9M M 

1 + 7A 2 * pA °-Jr 



+ 



R 5 



WR 3 

l + ^ P A + O(R- 2 )). 



(28a 2 - 14a 2 - 9a 2 ) + 0(i?" 



(3.36) 



Before proceeding further, we comment on the mass of the system and definition of the 
center of mass for each star. This is because there are several definitions of them in the PN 
approximation, and we should clarify the difference between the similar ones. 

First, we consider the conserved mass which is defined as 



M * = j d 3 xp 
= M 



c 2 V AR 
Using M*, f2 2 becomes 



1 /13M 12vrpAo M , n . 2 „ 2 , r 2N „ 
1 + — ( + - + ^rr^(34a 2 - 11a 2 , - 15a 2 ) + 0(iT 5 ) 



20R 3 



(3.37) 



tti™ = 3Ma?/35, Jiijj = Ma 2 a 2 /35 (i ^ j) 
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n 2 



2M 
R 3 



i + i 



11M, 
2# 



M* 
20.R 3 ' 



(90a? - 39a^ - 33a^) I + 0(i? 



18(^i 



i j* 



i? 5 



1 f 8 , 13M* ^„_ 2 , 



AR 



(3.38) 



where ( ^n)* = M*-t\\/M. Thus, fl 2 looks as if it depends on the internal structure of the 
star even for the limit di/R — > 0. Since we believe that in the equation of motion(EOM) for 
the point particle, the quantities depending on the internal structure does not appear, M* 
is not desirable to describe the EOM for the point particle. Instead, in the case when the 
EOM is derived, one usually adopts the PPN mass defined as 



Mi 



PPN 



/ d 3 xp 









V 

~2 



; self 



M 



1 /13M , M , 9 „ ,. „ 

1 + — I — ttt" + 2ixpA Q + 7^7^(38a? - 7aj - lha\) + 0{R °] 



(3.39) 



V AR ' ~" r " u ' 20.R 3 ' 

where U se \f and u se if are the self gravity part of the Newtonian potential and the spin velocity 
of each star. When we rewrite Eq.(|3.36|) using the PPN mass, the orbital angular velocity 



does not depend on the internal structure of the star and agrees with that for the point 
particle Q] in the limit a^/R — > as 



n 2 



2Mi 



PPN 



1 

1 + - 

cr 



R 3 

18(^~ii)ppn 



UMi 



PPN 



Mi 



PPN 



2R 



20R 3 



(94a? - 35a^ - 33a^) + 0(R~ 



R 5 



1 + 



-npA 



13M 



PPN 



+ 0(R~ 



(3.40) 



15 r AR 

where (^ii)ppn = Mppn^-ii/M. Thus, when we compare the present results with the point 
particle calculations, we should use the PPN mass. In the present case, however, Mppn is 
not a conserved quantity although M* is. When we consider a sequence of the equilibrium 
configuration as an evolutionary sequence, we should fix M*. 

Next, we consider the definition of the center of mass for each star. In the PPN formalism, 



it is defined as |24| 



x 



PPN 



M 



PPN 



d 3 xpx l 



l + 4f- + 3C/-V self + ^ 
c 2 2 2 2 



(3.41) 



and the x\ coordinate of the center of mass for star 1 deviates from to 
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Thus, in the PPN formalism, the following orbital separation should be used; 

AMaj 



R 



PPN 



R 



1 

1 + ^ 



5R 3 



+ 0(R~ 



(3.43) 



It is worth noting that when we define the center of mass by the conserved mass as 

d?xp«x\ (3.44) 



1 



the result is the same up to 0(i?~ 4 ). Thus, in this paper, we do not have to distinguish i?* 
from -Rppn- Even in the general cases, the difference between R* and -Rppn is expected to 
be small. 

Using R* and/or -Rppn, ^ 2 is rewritten as 



fi 2 



2M t 



l | 1 f 27rpA 11M* 



18(^u), 
R 5 , 



1 



1 + - \ -tftpA 



2R* 20R3 
13M* 



^(lSSa^-Sg^-SSa^j+O^; 4 ) 



4iL 



+ 0(R~ 



(3.45) 



or 



2M, 



PPN 



R3 

-flpp N 



18(^i 



1 + 



1 J PPN 



UMi 



PPN 



Mi 



PPN 



2R 



PPN 



20i?pp^j 



142a 2 - 35a 2 - 33a 2 ) + 0(R^ N ) 



R 5 



PPN 



, 1 f 2 „ 



13M P 



+ 0(i? 



PPN/ 



(3.46) 



4-Rppn 

Here, we should note that the effect of the spin-orbit coupling terms will appear in Q 2 from 



0(-R P p N ) |§. According to the PN study, it becomes 



a 2 



2M 



PPN 



C>3 

-flpp N 



1 - 



2M 



PPN / 2 



2D3 



<h + a 2 ) 



PPN 



(3.47) 



where we omit other terms which do not concern in this discussion. Eq.( p.47|) shows that 
the terms of 0(_Rpp N ) in Eq. (|3.46 ) cannot be explained only by the spin-orbit coupling term. 
This means that there is a new effect, say the PN quadrupole one, in Eq. (|3.46|) . 



19 



IV. THE ENERGY AND THE ANGULAR MOMENTUM 



A. The Total Energy 



The PN total energy is calculated from E = E^ ei + E^/c 2 , where g |0 



For the PN Darwin problem, they become 

12 D2n2 



E 



del 



M 



— + yK + i 2 ) + — — + °(^ ) 



= 2M 



l -^pf{llAl + a\A\ + 4Al + aiAl) 



R 3 
AMnpAo 
R 



Mnp / 87 4 4 4 \ 
5 [M 2 



M 2 

2 I R' 2 + ^( 5a l ~ 2a l ~ 2a D 



+tt 2 I— MR + (6a? + 21a 2 - 17 a: 

1 16 80R K 1 2 

{R 2 3 
— {3A - a\A 2 ) + -A {al + a 2 ) 



12 

35 1 



a\A x + a 2 A 2 



-—aja 2 2 (A 1 + A 2 )j 

+ l±n 2 (— + -(a 2 + a 2 2 )\ + -nU — + —(3a 2 + a 2 U + 0(iT 
p [ 10 35 V 1 2; J 8 1 16 10 V 1 - ' 



Pn 



i? 4 P 2 



2M 



34 / a ^2 

- — (npAo) 2 - 



llMnpA 7M 2 MnpA f 68 



3R 



32R 2 



+ 



R 3 



105 



(a 2 + a 2 ) 



61 



ii 



7 M 



M 2 



240P 4 



(302a 2 + 59a 2 - 209a 2 ;) + 0(R~ 



(4.1) 
(4.2) 

(4.3) 



(4.4) 



(4.5) 



where we use the equilibrium equations ( f2.12|) — ( |2.15| ) and Q 2 to reduce Eq. (|4.4|) to Eq.( |4.5| ). 

If we substitute Q 2 into the above formulas, E may be rewritten as -En + -^pn / c 2 where 
AnpAo 



£ N = M 
E PN = 2M 



5 
34 
21 



(TTpAoY 

M 2 
+ 240P 4 



M _J_ 3j " n i / 2 , n 2\ , n(R-5\ 

19M7rpA 25M 2 M7rpA f 22 



(4.6) 



6i? 



32i? 2 



+ 



i? 3 



2 2 158 ^ n 
2r ai + a2) -^^ 



(26a 2 + 35a 2 - 155a 2 ) + 0{R 



-5\ 



(4.7) 
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B. The Total Angular Momentum 



We can calculate the PN total angular momentum from J = + Jp^f/c 2 , where [22 



j: 



clef 



def 
PN 



pVtpd 3 X 
P 



v v iv z + 6U + - 



and Vu, = flzu 2 . For the PN Darwin problem, they become 



j; 



R 2 aj + af 



= MVt 



R 2 npA + 



Po \R 2 



4 I R 2 



p I 5 35 v 1 z/ | 5 



47rp 
1}5~ 



{l8A (a 2 + a 2 ) - YZ(a\A x + a\A 2 ) - a^a 2 ^ + A 2 )} 

M 



+ T + T { 1 + a ') + — r~ 



+ 



20R 



10a 2 + 27a 2 - 19a 2 ) 



MO* 



164 

i? 2 7rpA + 5RM H vrpA (a 2 + a 

105 



M 



10i2 



(20a 2 + 15a 2 , - 11a 2 ) + 0{R' 



(4.8) 
(4.9) 

(4.10) 



(4.11) 



(4.12) 



where we use the equilibrium equations (|2.12|) — (|2.15|) and Q 2 to reduce Eq. ( |4.11| ) to 

If we substitute Q 2 into the above formulas, J may be rewritten as Jn + Jpn / c 2 where 

R 2 a 2 + a 2 , \ 



2Mfi N h 



Ji 



PN 



v 4 

^ 2 7rpA + ^i?M + ^(2018a 2 + 2081a 2 + 21a 2 ) 



M 
80R 



122a 2 + 85a 2 . - 97a 2 ) + 0{R- 



(4.13) 



(4.14) 



V. EQUILIBRIUM SEQUENCE OF THE POST-NEWTONIAN DARWIN 

ELLIPSOID 

In this section, we construct the equilibrium sequences of the PN Darwin ellipsoids fixing 
M* and p as follows; 
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(1) using Eqs.( p.l3|) — ( |2.15|) , we numerically calculate the equilibrium sequences of the New- 
tonian order. Up to this stage, a%, a 3 , and R = R/ai are determined. 

(2) di is determined from the condition M* =constant using Eq. fl3.37|) : 



3M* 



47T pa 2 a 3 



l _ up/ 3M> \ 2 / 3 /12Ap 13a 2 a 3 a 2 q 3 ^ 
c 2 Uvrpa 2 a 3 / V 5 3.R 15.R 3 



2 - 15o£ 



• (5-1) 



(3) after substituting Eq.flO]) into Eqs.(CT), flO|) , (gg ), (^TD| ), and we rewrite the 

PN expressions for the orbital angular velocity, the energy and the angular momentum as 



n 2 

E 
J 



^ - f2 2 + ^ 2 



(M 2 /a*) 
J 

(MX) 1/2 



-En H — ^ — -Epn, 

- Jn H 5 ^PN, 

era* 



where 



a* = 


/3M*\ 


1/3 


2 - 


4 


■ 2 

-I 3 


o 2 - 


^("2«3 


N4/3 



3 A 



9 42 A , 



R 3 a 2 a 3 2i? 4 25R 5 a 2 a 3 



1 



(a 2 a 3 ) 



1/3 



5i? 6 
3 Ac 



^(28-14a^-9a^ 



1 



5 a 2 a 3 2R 10R 3 



[2 - ad 



3f2^j 



20a 2 a 3 



E PN = (a 2 a 3 



.2/3 



Ac 



140 \a 2 a 3 
1 



55 

+ =T + 



1 Ar 



JpN 



2(a 2 a 3 ) 



-1/6 



120i? 4 

3fi^j 
4a 2 a 3 



48i? 2 700.R 3 a 2 a 3 
194 + 215a^ + 165a^) 



(398 + 401a 2 + a 



1/2 



4 



(a 2 a 3 ) 



+ 



1/6 



An 



3^^ 


1/2 . 


_4a 2 a 3 _ 





8 a 2 a 3 48 



1400 a 2 a 3 



(338 + 401a 2 + 21a 2 



240R 



-(494 + 155a 2 + 141a 2 ) 



(5.2) 
(5.3) 
(5.4) 

(5.5) 
(5.6) 



(5.7) 



(5.8) 



(5.9) 
(5.10) 



(5.11) 
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Then, using the numerical value of a>2, «3, and R determined at (1), we calculate the sequence 
of the angular velocity, the energy and the angular momentum as a function of the orbital 
separation. 

We repeat this procedure changing the mean radius of each star a*. Once a sequence is 
obtained, we search the minimum point of the energy. If we find it, we call it the ISCCO. 

In figs. 2(a) and (b), we show E and J as functions of the normalized separation R*/a*, 
where E = E/(M*/a,) and J = J/{M*a*) 1 / 2 . The fi gures show the important fact that the 
orbital separation at the ISCCO decreases approximately in proportion to M*/c 2 a*, i.e., the 
characteristic value of the compactness of each star. 

In fig. 3, we show the normalized orbital angular velocity Q = Q/^yWp at the ISCCO as 
a function of compactness. We show Q at the energy minimum as well as that at the angular 
momentum minimum. The figure indicates that two minimums are almost coincident, but 
slightly different. The deviation between the location of the minimums of the energy and 
the angular momentum comes from the fact that we assume a*/i?* is a small parameter, and 
expand the energy to 0(R~ A ) and the angular momentum to 0(R~ 1 ). In any case, we may 
expect that the ISCCO locates near two minimums. Fig. 3 clearly shows that the orbital 
angular velocity at the ISCCO increases almost linearly with increase of the compactness. 
Since the PN approximation is valid only for small compactness, we can not do solid estimate 
of Q at the ISCCO for a realistic NS of M*/c 2 a* ~ 0.2. Here, we dare to extrapolate this 
results to the realistic NS's. Then, we can find that the angular velocity at the ISCCO is 
about 10% larger than that in the Newtonian case. 

Before closing this section, we summarize the numerical results in Table I. In Table I, 
f denotes the point of the ISCCO defined by the energy minimum. We note that in the 
Newtonian order, the axial ratios «2 and at the ISCCO are 0.8573 and 0.7995, respectively 
|23|j . Thus, the PN effect increases c*2 and 0:3 at the ISCCO; i.e., tidal deformation is weak 
compared with the Newtonian case. The reason for such a behavior seems to be that each 
star of binary is forced to be compact due to the PN gravity (see Eq.( |5.1| )), and as a result, 
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the tidal force is less effective than in the Newtonian case. 



VI. SUMMARY 

In this paper, we have calculated the equilibrium sequences of the corotating BNS's of the 
incompressible fluid using the first PN approximation of general relativity. The conclusions 
are as follows. 

(1) due to the PN effect, the orbital separation at the ISCCO(secular instability limit) 
decreases in proportion to the compactness of the star M^/c^a*. 

(2) the orbital angular velocity at the ISCCO increases in proportion to the compactness of 
the star M*/c 2 a*. 

(3) the reason for features (1) and (2) is that each star is forced to be compact due to the PN 
gravity and as a result, the tidal effect becomes less effective compared with the Newtonian 
binary. 

These results agree with the recent numerical study by Shibata [^(J . 

Since the analysis in this paper is done almost analytically, the equilibrium sequence ob- 
tained may be regarded as the exact solution of the Einstein equation for the limit of small 
compactness M*/c 2 a* <C 1 and a* <C R. Recently, two groups have been performed numer- 
ical computations for obtaining the equilibrium state of BNS's using the semi-relativistic 



approximation [25]. The present solution is useful to check their numerical solutions as 
well as to understand the essence in their results. 

Finally, we comment on the possibility of extending this work to the compressible fluid 
case. In this paper, it is found that the deformation of the ellipsoidal figure due to the PN 
gravity is small effect to the angular velocity, so that the ellipsoidal model for the density 
configuration may be a good approximation. If we determine the density configuration as 
the ellipsoidal approximation as was done by Lai, Rasio and Shapiro ||23|| , we only need to 
carry out the integrals which appear (a) in calculating the angular velocity by the first TV 
equation , and (b) in calculating the energy and the angular momentum, which can be done 



24 



easily. We will publish results of such a calculation in a subsequent paper 
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TABLE CAPTIONS 

Tabled. Equilibrium sequences of the PN Darwin ellipsoids with M^/c 2 a if = 0.01, 0.03 
and 0.05. f denotes the point of the ISCCO. 
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TABLES 



R*/a* 


#*/(«!)* 


a 2 /ai 


a 3 /ai 


n 


J 


E 


M*/c 2 a* = 0.01 


2.80 


2.190 


0.7283 


0.6571 


0.3645 


1.574 


-1.303 


3.00 


2.542 


0.8135 


0.7483 


0.3208 


1.527 


-1.317 


3.202f 


2.829 


0.8593 


0.8020 


0.2877 


1.518 


-1.320 


3.25 


2.891 


0.8673 


0.8119 


0.2809 


1.518 


-1.320 


3.50 


3.205 


0.8999 


0.8532 


0.2497 


1.528 


-1.317 


4.00 


3.785 


0.9372 


0.9040 


0.2033 


1.571 


-1.309 


M*/c 2 a* = 0.03 


2.80 


2.296 


0.7776 


0.7087 


0.3488 


1.534 


-1.311 


3.00 


2.602 


0.8388 


0.7775 


0.3101 


1.510 


-1.318 


3.123| 


2.768 


0.8633 


0.8069 


0.2906 


1.507 


-1.319 


3.25 


2.931 


0.8828 


0.8312 


0.2729 


1.509 


-1.318 


3.50 


3.235 


0.9106 


0.8673 


0.2433 


1.521 


-1.315 


4.00 


3.805 


0.9433 


0.9126 


0.1987 


1.565 


-1.307 


M*/c 2 a* = 0.05 


2.80 


2.370 


0.8121 


0.7467 


0.3353 


1.508 


-1.314 


3.00 


2.651 


0.8596 


0.8024 


0.2999 


1.496 


-1.318 


3.044f 


2.708 


0.8674 


0.8119 


0.2932 


1.495 


-1.318 


3.25 


2.966 


0.8962 


0.8484 


0.2649 


1.500 


-1.316 


3.50 


3.262 


0.9201 


0.8801 


0.2368 


1.514 


-1.313 


4.00 


3.824 


0.9488 


0.9206 


0.1939 


1.558 


-1.305 



TABLE I. 
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FIGURE CAPTIONS 



Fig.l. Sketch of the Darwin ellipsoid. The origin of the coordinate we choose locates at 
the center of mass of star 1. 

Fig. 2(a). The total energy of the equilibrium sequence as a function of R*/a*. Solid, 
dotted, dashed and long-dashed lines denote M*/c 2 a* = 0(the Newtonian case), 0.01, 0.03 
and 0.05, respectively. 

Fig.2(b). The total angular momentum of the equilibrium sequence as a function of 
R*/a*. Solid, dotted, dashed and long-dashed lines denote M*/c 2 a* = 0(the Newtonian 
case), 0.01, 0.03 and 0.05, respectively. 

Fig. 3. The orbital angular velocity at the ISCCO as a function of the compactness 
parameter M*/ c 2 a*. Filled circles and open triangles denote fl/^/Wp for the minimum points 
of the total energy and total angular momentum, respectively. 
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